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We show for the Ising model that is possible construct a discrete time stochastic model analogous 
to the Langevin equation that incorporates an arbitrary amount of damping. It is shown to give 
the correct equilibrium statistics and is then used to investigate nonequilibrium phenomena, in 
particular, magnetic avalanches. The value of damping can greatly alter the shape of hysteresis 
loops, and for small damping and high disorder, the morphology of large avalanches can be drastically 
effected. Small damping also alters the size distribution of avalanches at criticality. 
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I. INTRODUCTION 

In many situations, it is useful to discretize continu- 
ous degrees of freedom to better understand them, both 
from a theoretical standpoint and for numerical effi- 
ciency. Ising models are perhaps the best example of this 
and have been the subject of numerous theoretical and 
numerical studies. Renormalization group arguments 1 
have explained the reason why this discretization gives 
equilibrium critical properties of many experimental sys- 
tems, and these kinds of arguments have been extended 
to understanding their equilibrium dynamics 2 . For non- 
equilibrium situations, such as the study of avalanches, 
such arguments probably do also apply to large enough 
length and time scales as well. However there are many 
situations where it would be desirable to understand 
smaller length scales where other factors should become 
relevant. 

This is particularly true with dynamics of magnetic 
systems, where damping is often weak in comparison to 
precessional effects. For studies of smaller scales, it has 
been necessary to use more time consuming micromag- 
netic simulations utilizing continuous degrees of freedom, 
such as the Landau Lifshitz Gilbert equations 3 which is a 
kind of Langevin equation that gives the stochastic evo- 
lution of Heisenberg spins. 



ds 

dt 



-s x (B — 7s x B), 



(1) 



where s is a microscopic magnetic moment, B is the local 
effective field, and 7 is a damping factor, measuring the 
relative importance of damping to precession. In real 
materials it ranges- from small damping 7 = .01, to 1. 
In contrast, the dynamical rules implemented for Ising 
models are most often "relaxational" so that energy is 
instantaneously dissipated when a spin flips, as with the 
Metropolis algorithm. 

However there is a class of "microcanonical" Ising dy- 
namics^ reviewed in section UT1 where auxiliary degrees of 
freedom are introduced and all moves conserve the total 
energy. The other degrees of freedom can be taken to be 



variables associated with each spin, and allowed moves 
can change both the state of the spins and the auxiliary 
variables. This can be thought of crudely, as a discretized 
analogy to molecular dynamics, and is also similar to dis- 
crete lattice gas models of fluids 6 -' 7 ". These models give 
the correct equilibrium Ising statistics of large systems 
and can also be used to understand dynamics in a differ- 
ent limit than the relaxational case. 

Real spin systems are intermediate between these two 
kinds of dynamics and as mentioned above, are better de- 
scribed by Langevin dynamics. In the context of spins, 
the question posed and answered here is: how does one 
formulate a discrete time version of stochastic dynamics 
that includes damping and gives the correct equilibrium 
statistics? In section Hi Al we are able to show that there 
is a fairly simple method for doing this using a combina- 
tion of microcanonical dynamics, and an elegant proce- 
dure that incorporates damping and thermal noise. This 
procedure differs from that of the Langevin equation in 
that it requires non-Gaussian noise. Despite this, the 
noise has surprisingly simple but unusual statistics. 

We will then show that this procedure gives the correct 
equilibrium statistics and verify this numerical in section 
IIIBI with a simulation of the two dimensional Ising model 
with different amounts of damping. 

Because the value of damping is an important phys- 
ical parameter in many situations it is important that 
there is a straightforward way of incorporating its effects 
in Ising simulations. This is particularly noteworthy as 
Ising kinetics are a frequently used means of understand- 
ing dynamics in many condensed matter systems. 

After this in section HTT1 we will turn to nonequilibrium 
problems where, using this approach, we can study the 
effects of damping on a number of interesting properties 
of systems displaying avalanches and Barkhausen noised. 
We first show how to modify the kinetics for this case and 
then study systems in two and three dimensions. With 
modest amounts of computer time, we can analyze prob- 
lems that are out of the reach of micromagnetic simula- 
tions and allow us to probe the effects of damping on the 
properties of avalanches. This is related to recent work 9 
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by the present authors using both the Landau Lifshitz 
Gilbert equation, Eqn. [TJ and theoretical approaches, to 
understand how relaxational dynamics of avalanches 10 , 
are modified at small to intermediate scales by this more 
realistic approach. With the present approach we find 
new features and modifications of avalanche dynamics. 
We find that the shape of hysteresis loops can be strongly 
influenced by the amount of damping. One of the most 
striking findings is that there exists a parameter regime 
of high disorder and small damping where single system- 
size avalanches occur that are made up of a large number 
of disconnected pieces. We can also analyze the criti- 
cal properties of avalanches when damping is small and 
give evidence that there is a crossover length scale, below 
which avalanches have different critical properties. 

II. NON-RELAXATIONAL DYNAMICS 

We start by considering a model for a magnet with con- 
tinuous degrees of freedom, such as a Heisenberg model 
with anisotropy. The Ising approximation simplifies the 
state of each spin to either up or down, that is Sj = ±1, 
i = 1, . . . , N. One important effect that is ignored by 
this approximation is that of spin waves that allow the 
transfer of energy between neighbors, and for small os- 
cillations, give an energy contribution per spin equal to 
the temperature T (here we set ks = !)• This moti- 
vates the idea that there are extra degrees of freedom 
associated with every spin that can carry (a positive) 
energy ej. Creutz introduced such degrees of freedom 5 
and posited that they could take any number of dis- 
crete values. He used these auxiliary variables to 
construct a cellular automota to give the correct equi- 
librium statistics for the Ising model, in a very efficient 
way that did not require the generation of random num- 
bers. Thus we have a Hamiltonian H to t that is the sum 
of both spin H spin and auxiliary degrees of freedom H e : 
Htot = H spin + H e . H spin can be a general Ising spin 
Hamiltonian and H e = X^ e «- in our model there is a 
single auxiliary variable e, associated with each lattice 
site i, that can take on any real value > 0. 

However for the purposes of trying to model dynamics 
of spins, it also makes sense to allow the e^'s to interact 
and exchange energy between neighbors. For example, 
one precessing spin should excite motion in its neighbors. 
This exchange was formulated in the context of solidifi- 
cation using a Potts model instead of an Ising model by 
Conti et aZ<ii, but can equally well be used here. 

Now we can formulate a microcanonical algorithm for 
the Ising model using a procedure very similar to their 
prescription. In each step: 

1. We choose a site i at random. 

2. We randomly pick with equal probability either a 
spin or an auxiliary degree of freedom, Sj of e%: 



(a) s^s: We attempt to move spins (such as the 
flipping of a single spin) . If the energy cost in 
doing this is < ej we perform the move and 
decrease e, accordingly. Otherwise we reject 
the move. 

(b) ej's: We pick a nearest neighbor j, and repar- 
tition the total energy with uniform probabil- 
ity between these two variables. That is, af- 
ter repartitioning, e[ = (e^ + ej)r and e'j = 
(e-i + ej)(l — r), where < r < 1 is uniform 
random variable. 

Note that these rules preserve the total energy and the 
transitions between any two states have the same proba- 
bility. Therefore this will give the correct microcanonical 
distribution. For large N, this is, for most purposes 12 , 
equivalent to the canonical distribution oc exp (—j3Htot)- 
Note that the probability distribution for each variable 
ej, P(e,) = /3exp (— /3ei), so that the (ej = T. That is, 
measurement of average of e^'s directly gives the effective 
temperature of the system. 



A. Extension To Damping 

The question we asked, is how to extend this equilib- 
rium simulation method to include damping. In this case 
the system is no longer closed and energy is exchanged 
with an outside heat bath through interaction with the 
auxiliary variables. As with the Langevin equation, there 
are two effects. The first is that the energy is damped. 
Call the dissipation parameter for each step a, which will 
lie between and 1. Then at each time step we lower 
the energy with e, — > aei for all sites i. By itself, this 
clearly will not give a system at finite temperature and 
we must also include the second effect of a heat bath, 
which adds energy randomly to the system. In the case 
of the Langevin equation, a Gaussian noise term n(t) is 
added to keep the system at finite temperature. A dis- 
cretized version of this, that evolves the energy e(t) at 
time step t is 

e(t+l) = ae(t) +n{t). (2) 

This equation will not work if the noise n(t) is Gaussian 
as this does not give the Gibbs distribution P eq (e) = 
(3 exp (— 0e). Therefore we need to modify the statistics 
of n(t). It is possible to do so if we choose n(t) at each 
time t from a distribution 

p(n) = aS(n) + (1 - eO/3e -/5n 0(n) (3) 

where 9 is the Heaviside step function. To show this, we 
write down the corresponding equation for the evolution 
of the probability distribution for e: 
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P(e',t + 1) = {5(e' - (ae + n))) = / P(e,t)p(n)5(e' - (ae + n))de dn 



(4) 



r 



We require that the P as t - 
P(e,t) =P eq {e) =/?exp(-/3e). 



* oo obeys P(e,t + 1) = 
for e > 0. It is easily ver- 
ified that by choosing this form of P(e, t) and by choosing 
P{n) as in Eq. [3l we satisfy Eq. [4] 

Therefore to add damping to this model, we add the 
following procedure to the steps stated above: 

3. Choose a uniform random number < r < 1. If r < 
a, then — > aej. Otherwise — > ae^ — Tln(r'), 
where r' is another uniform random number be- 
tween and 1. 

If we assume that the probability distribution for the 
total system is of the form Poibbs oc exp (—f3H to t) — 
cxp (— (3 H S pi n ) exp (— (3H e ), we will now show that the 
steps 1, 2, and 3, of this algorithm preserve this distri- 
bution. Following the same reasoning as above for the 
microcanonical simulation, moves implementing steps 1 
and 2 do not change the total energy, and they preserve 
the form of PGibbs because PGibbs depends only on the 
total energy (H to t), and 1 and 2 explore each state in 
an energy shell with uniform probability. Because of 
the form of Pcibbs, its dependence on the variable is 
ex exp (— I3ei). According to the above argument, after 
step 3, it will remain unchanged. Therefore all steps in 
this algorithm leave Poibbs unchanged. The algorithm is 
also ergodic, and therefore this will converge to the Gibbs 
distribution-^ as t — ► oo. 

Because the steps each preserve the Gibbs distribution, 
the ordering of them is not important in preserving equi- 
librium statistics. For example, we could sweep through 
the lattice sequentially instead of picking i at random. 
We could perform step 3 after steps 1 and 2 were per- 
formed N times. 



a = 0.9 but are so close as to be indistinguishable and 
are therefore not shown. We also checked that the distri- 
bution of auxiliary variables had the correct form. The 
probability distribution for the energy e is shown in Fig. 
[3] Fig. [5] plots the distribution P(ej) versus energy e», 
averaged over all sites i on a linear-log scale for a = 0.5 
and T = 0.8, and 1.1. The curves are straight lines over 
four decades and show the correct slopes, for T = 0.8, 
(e») = 0.8002 and for T = 1.1, (e,) = 1.1003. 
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FIG. 1: Plot of results obtained for the two dimensional Ising 
model on a 128 2 lattice for two different values of the damping 
parameter. This is a plot of the average magnetization per 
spin m vs. T. The x's are for no dissipation, a = 1, which is 
a purely microcanonical simulation. The +'s are for a = 0.5. 
The dashed curve is the exact solution to this model in the 
thermodynamic limit. 



B. Equilibrium Tests 

We performed tests on this algorithm and verified that 
it did indeed work as expected. We simulated the two 
dimensional Ising model on a 128 2 lattice with differ- 
ent values of the damping parameter, and compared it 
with the exact results. The average magnetization per 
spin m is plotted in Fig. [1] as a function of the tempera- 
ture T and compared with the exact result" for large N 
(dashed curve). The x 's are the case a = 1, which is then 
just an implementation of the microcanonical method^ 
described above. In this case, the temperature was ob- 
tained by measuring (a) because the energy was fixed 
at the start of the simulation. The only point which is 
slightly off the exact solution is in the critical region, as 
is to be expected. The case a = 0.5 is shown with the 
+'s and lie on the same curve. Results were obtained for 



III. AVALANCHE DYNAMICS 

Avalanche dynamics of spin systems have been mainly 
studied using models that are purely relaxational. There 
is a whole range of interesting phenomena that have been 
elucidated by such studies and have yielded very inter- 
esting properties. The simplest model that can be used 
in this context is the random field Ising model (RFIM) 
with a Hamiltonian 



H 



jsjSj - kiSi -hy~]sj 



(5) 



where J is the strength of the nearest neighbor coupling, 
hi is a random field, with zero mean, and h is an exter- 
nally applied field. A magnet is placed in a high field h 
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FIG. 2: Plot of results obtained for the two dimensional 
Ising model on a 128 2 lattice for the probability distribution 
for the auxiliary variables a, at two different temperatures 
with a damping parameter a = 0.5. The upper curve is for 
T = 1.1 and the lower for T = 0.8. 

and then this is very slowly lowered. As this happens, the 
spins will adjust to the new field by flipping to lower their 
energy. In the usual situation, the system is taken to be 
at T = 0, so that only moves that lower the energy are 
accepted. The flipping of one spin can cause a cascade of 
additional spins to flip, causing the total magnetization 
M to further decrease. The occurrence of these cascades 
is called an "avalanche" . At zero temperature there is 
one parameter j that characterizes the system, the ratio 
of nearest neighbor coupling to the distribution width of 
the random field. One considers the behavior of a sys- 
tem when its starts in a high field and is slowly lowered. 
When j is small the system is strongly pinned and the 
system will have a number of small avalanches generat- 
ing a smooth hysteresis loop. For large j, the system 
will have a system-size avalanche involving most of the 
spins in the system, leading to a precipitous drop in the 
hysteresis loop. There is a critical value of j where the 
distribution of avalanche sizes is a power law and self- 
similar scaling behavior is observed. 

Here we investigate how this is modified by adding 
damping to these zero temperature dynamics according 
to the following rules: 

1 . The field is slowly lowered by finding the next field 
where a spin can flip. 

2. The spins then flip, exchanging energy with auxil- 
iary variables e, as described above. The number of 
times this is attempted is n m times the total num- 
ber of spins in the system. Here we set n m = 16. 
In more detail: 

(i) Spin moves: An attempt to move each spin 
on the lattice is performed by attempting to 
flip sequentially every third spin, in order to 



minimize artifacts in the dynamics due to up- 
dating contiguous spins. (The lattice sites are 
linearly ordered using "skew" boundary con- 
ditions). Then all three sublattices are cycled 
over. 

(ii) Energy moves: Exchange of energy 
with nearest neighbors is performed cycling 
through all directions of nearest neighbors. 
Using the same sequence of updates, the e^'s 
exchange energy with their nearest neighbors 
in one particular direction. 

(iii) Dissipation: The energy of each ej is lowered 
to aej. 

3. We check for when the spins have settled down as 
follows: if the e.;'s are not all below some energy 
threshold ethresh, set below to be 10 -4 , or the spin 
configuration has changed, step 2 is repeated until 
these conditions are both met. 

4. When the spins have settled down, we go to step 1. 

The parameters n m and ethresh were varied to check that 
the correct dynamics were obtained. The larger a, the 
smaller the dissipation and the larger the number of iter- 
ations necessary to achieve the final static configuration. 
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FIG. 3: The major branch of the descending hysteresis loop 
for 64 2 systems using different values of the damping parame- 
ter and the spin coupling. Strong damping, a = 0.5 is shown 
in the left most curve (as judged from the top of the plot) 
for coupling j = 0.3 which starts decreasing from M — 1 at 
h — —0.2, and does not have large abrupt changes. All the 
other curves are for weak damping, a — 0.99. In this case but 
also for j — 0.3, we see that although M starts to decrease at 
the same location as for strong damping, it drops abruptly as 
the field is lowered. As the coupling j is decreased, smooth 
curves are eventually seen again. Going left to right, as judged 
from the top, are j = 0.3, 0.25, 0.2, and 0.15. 



(a) 



(b) 



FIG. 4: (a) The spin configuration for a 256 2 system with j = 0.35, a — 0.9 during a system size avalanche at the field 
h — —0.400007. (b) A gray-scale plot of the auxiliary variables at the same time. 




FIG. 5: Spin configurations for a 256 2 system with j = 0.25, a — 0.99 during an avalanche at the field h = — 7 x 10~ 5 . (a) 
The beginning of the avalanche, (b) When the avalanche is of order of half the system size, (c) The final configuration of the 
avalanche. 



A. Two Dimensional Patterns 

We first investigate the case of two dimensions where 
it is simpler to visualize the avalanches in various condi- 
tions than in three dimensions. Much experimental work 
and theoretical work on avalanches has been done on two 
dimensional magnetic films and this case should be highly 
relevant^. 

We first examine how the hysteresis loops change as 
a function of the coupling j and the damping parame- 
ter a for a 64 2 system. The major downwards hystere- 
sis loops are shown in Fig. [3] for a variety of param- 
eters described below. We first examine strong damp- 
ing a — 0.5. For j — 0.3 the hysteresis curve is quite 
smooth with all avalanches much less than the system 



size (left most curve) . Now consider the same value of j 
but with with small damping, a — 0.99. The curve now 
is a single downwards step with a small tail at negative 
h. The lower damping has allowed that system to form 
a system size avalanche. The difference is due to the 
fact that with small damping, the energy of avalanched 
spins is not immediately dissipated and as a consequence, 
heats up neighboring spins, allowing them to more easily 
avalanche as well. Therefore a system size avalanche is 
seen in the small damping case, leading to the precipitous 
drop in the hysteresis loop. 

When the value of the coupling j is lowered to 0.15 for 
a = 0.99, smooth loops are obtained. The Fig. [3] shows 
intermediate values of the coupling parameter as well. 

To better understand the reason why the energy of 
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the auxiliary variables can trigger further spins to flip, in 
Fig. Q]we show the state of a system during a system size 
avalanche for j = 0.35 and a moderately small damping 
value, a = 0.9, with h = -0.400007. Fig. gfa) shows that 
the flipped spins form a fairly compact cluster and Fig. 
IDJb) shows the corresponding values of the e,'s in a gray 
scale plot, suitably normalized. It has the appearance of 
a halo around the growth front of the avalanche. The 
spins in the growth front have just flipped and so energy 
there has not had a chance to diffuse or dissipate and 
so has a higher spin temperature. The interior is cold 
because damping has removed energy from the auxiliary 
degrees of freedom. This higher temperature diffuses into 
the the unflipped region allowing spins to flip by thermal 
activation. 

Because large avalanches are possible for small damp- 
ing in a parameter range where the relative effect of the 
random field is much larger, it is of interest to see if 
avalanches have a different morphology than typical large 
avalanches for high damping systems. Fig. shows such 
spin configurations first at the beginning of the avalanche 
and further along during propagation when it has reached 
roughly half the system size, and finally when it has 
reached its final configuration and the maximum aux- 
iliary variable value is < 4 x 10~ 4 . The morphology of 
this is very different than what is seen for large avalanches 
with stronger coupling, for example Fig. [3J At very small 
fields, in this figure h = — 7 x 10~ 5 , surface tension pre- 
cludes the formation of minority domains, but because 
disorder is large, there will be many small regions where 
the local field is much stronger and these will want to 
form downward oriented (black) domains. There is a fi- 
nite activation barrier to forming these that can only be 
overcome at finite temperature. However the majority 
of the spins still strongly disfavor flipping. But because 
damping is small, heat has a chance to diffuse through 
these regions into the favorable regions, allowing discon- 
nected regions to change orientation by thermal activa- 
tion. Note that we have checked numerically that small 
damping with strong coupling also leads to compact con- 
figurations, so disorder is an essential ingredient in this 
new morphology. 



B. Three Dimensions 

We first check that as with two dimensions, the value 
of the damping parameter can have a large effect on the 
shape of a hysteresis loop. Fig. [6] shows the downward 
branches of the major hysteresis loop when the only pa- 
rameter that is changed is the damping, a. The system 
is a 32 3 lattice with j = 0.19. A value for high damping, 
a = 0.5, is the upper line. The lower line is for small 
damping with a — 0.99. 

A more subtle effect, is that of damping on what hap- 
pens near criticality. In this case the value of the criti- 
cal j will depend on the value of a as is apparent from 
the results of Fig. [5J At this point, the distribution of 




FIG. 6: The major branch of the descending hysteresis loops 
in two 32 3 systems with j = 0.19, for two different values of 
the dissipation, upper curve: a — 0.5, lower curve: a — 0.99. 



avalanche sizes is expected to follow a power law distribu- 
tion for large sizes. We located this point and examined 
system properties in this vicinity. Fig. [7] shows examples 
of such runs for 32 3 systems. Fig. [Tfa) shows a plot of 
the magnetization per spin M, versus the applied field h 
for j = 0.165 and j = 0.167. For larger values of j, the 
avalanches rapidly become much larger as is seen in Fig. 
IH1 and for smaller values, avalanches all become small. 
Fig. \T(b) shows a plot of the same quantity with relax- 
ational dynamics near criticality. The avalanches take 
place over a much smaller range in applied field. 

To quantify this difference, we studied the avalanche 
size distribution exponent that is obtained by calculating 
the distribution of avalanche sizes over the entire hys- 
teresis loop. This was studied by averaging avalanches of 
many runs, (200 for a — 0.99) for 32 3 systems and for dif- 
ferent values of parameters. We show a comparison of the 
avalanche size distribution for a — 0.99, shown with +'s 
and for a — 0.9, shown with x's in Fig. [5] For a = 0.99 
the curve fits quite well to a power law with an exponent 
of —1.4 ± .1 as shown in the figure. For purely relax- 
ational dynamics, the same exponent has been carefully 
measured^ to be 2.03 ± .03 (which is consistent with our 
results for relaxational dynamics on much smaller sys- 
tems than theirs). With smaller damping we expect to 
have a crossover length corresponding to the length scale 
associated with the damping time, above which the dy- 
namics should appear relaxational. a — 0.9 appears to 
show such a crossover from a slope of approximately —1.4 
for small avalanches, to a higher slope for large ones. A 
line with slope of —2 is shown for comparison and appears 
to be consistent with this interpretation. 
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(a) h (b) h 

FIG. 7: (a) Magnetization versus field for the Ising model with damping described in the text. The system size is 32 3 and 
the two lines represent two runs close to criticality, one with a coupling of j = 0.165 and the other of 0.167. (b) The plot for 
relaxational dynamics (large damping) with couplings of .21 and .212. 
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FIG. 8: The avalanche size distribution, measured of the 
entire hysteresis loop for a — 0.99 (+ symbols) and a — 0.9 (x 
symbols) . The x-axis is the number of avalanches normalized 
by it's mean size. The y-axis is the normalized distribution 
of sizes. The less negative sloped straight line is a fit of the 
a = 0.99 curve and has a slope of —1.4. The more strongly 
sloped one has a slope of —2. 



IV. DISCUSSION 

This paper has introduced a new set of dynamics for 
Ising models that incorporates damping in a way that 
has not before been achieved. The dynamics that have 
been devised have a lot in common with Langevin dy- 
namics, except they are for discrete rather than contin- 
uous systems. In Langevin equations, a continuous set 
of stochastic differential equations are used to model a 
system. It differs from molecular dynamics in that ther- 
mal noise and damping are both added so that the sys- 
tem obeys the correct equilibrium statistics. In the case 
studied here, we start by considering microcanonical dy- 



namics 5,11 which introduces auxiliary degrees of freedom. 
We then add damping and thermal noise. Whereas the 
thermal noise is typically Gaussian in the case of the 
Langevin equation, here it must be taken to be of a spe- 
cial exponential form, Eq. [3l in order for it to satisfy the 
correct equilibrium statistics. 

The form of this noise, although quite unusual, can 
be understood, to some extent qualitatively. For large 
damping, or small a, the strength of the 5 function be- 
comes small, and the effect is dominated by the second 
term which is oc exp(— /3n) (for positive n). Although this 
is non-Gaussian, n can be thought of as a random amount 
of positive energy. In the Langevin equation, noise is of- 
ten added to a velocity degree of freedom. In terms of 
a velocity, the exponential form that we have obtained 
would correspond to a Gaussian if this was expressed in 
terms of a velocity instead. In the limit of small damp- 
ing, where a is close to 1, the effect of the noise becomes 
small because the first term, which is to add no noise, will 
dominate the distribution. This is in accord with what 
happens in the Langevin equation where if dissipation is 
small, little thermal noise is needed to keep the system 
at a given temperature. 

The fact that it is possible to model damped systems 
in this discrete manner should have many useful applica- 
tions, and is easily extended to other kinds of systems, 
aside from Ising models, especially in applications where 
computational efficiency is an important criterion. 

The case of avalanches in magnetic systems is an in- 
teresting nonequilibrium use of these dynamics. Al- 
though one might expect that in most situations, for large 
enough distance and time scales, finite damping will be 
unimportant, physics at smaller scales is still of great in- 
terest and effects at those scales can propagate to larger 
scales. Because damping in real materials can be quite 
small, their effects are readily observable experimentally. 
This work is expected to be important at intermediate 
scales. We have investigated the phenomenon seen in this 
model with varying degrees of damping and found that 
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it makes a qualitative difference to many of the features 
seen on small and intermediate scales. This work is by no 
means exhaustive and there are many other effects that 
can be investigated by straightforward extensions. The 
effect of dipolar interactions in conjunction with damp- 
ing could also be explored. We have chosen to update the 
spin and auxiliary variables at equal frequencies. Varying 
this should lead to a different value for the heat diffusion 
coefficient which should change the quantitative values 
for length and time scales. 

The phenomena we have found was in qualitative 
agreement with earlier work using the Landau Lifshitz 
Gibbs equations^. As avalanches progress, the effective 
temperature, which we have seen can be quantified by 
(ej) at site i, will increase as energy is released. This 



energy then diffuses to the surrounding regions, giving 
those spins the opportunity to lower their energy by ther- 
mal activation. This allows avalanches to more easily 
progress when the damping is small in contrast to relax- 
ational dynamics, which has effectively infinite damping, 
a = 0. This can lead to some substantial differences in 
avalanche morphology, particularly as for small damp- 
ing, highly disordered systems can avalanche. At low 
fields this leads to a single avalanche being composed of 
many disconnected pieces. Experiments have been de- 
vised^ that are close experimental realization of the two 
dimensional random field Ising model, and it would be 
interesting to determine if systems such as this one, or 
similar to it, show avalanches with this morphology. 
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